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A Priori Error Estimates for Some Discontinuous Galerkin 
Immersed Einite Element Methods*^ 

Tao Lin J Qing Yang§ Xu Zhang^ 


Abstract 

In this paper, we derive a priori error estimates for a class of interior penalty discontinuous Galerkin 
(DG) methods using immersed finite element (IFE) functions for a classic second-order elliptic interface 
problem. The error estimation shows that these methods can converge optimally in a mesh-dependent 
energy norm. The combination of IFEs and DG formulation in these methods allows local mesh re¬ 
finement in the Gartesian mesh structure for interface problems. Numerical results are provided to 
demonstrate the convergence and local mesh refinement features of these DG-IFE methods. 


1 Introduction 


Let be a rectangular domain in and let F C D be a smooth curve separating fl into two sub- 
domains 0.~ and with n = 0 (see the first plot in Figure]^. We consider the following typical 
elliptic interface problem 

-V ■ iPVu) = f, in D+UD", (El) 

M = 0, on do, (1.2) 


where the diffusion coefficient /3 is a positive piecewise constant function: 

(‘- 3 ) 

According to the conservation laws, the following jump conditions are required on the interface: 


Mr 


/3 



0 , 

0 . 


(1.4) 

(1.5) 


Interface problems arise in many applications where mathematical simulations are carried out in a 
domain containing multiple materials. The elliptic interface problem (1.11 - (1.51 considered in this article 
appears frequently because the involved differential equation captures many basic physical phenomenons. 
A wide variety of numerical methods have been developed for interface problems, among which the finite 
element methods are advantageous for their capability to handle simulation domains with complicated 
geometry. It is well-known that conventional finite element methods generally require the mesh to fit the 
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interface geometry (see the second plot in Figure [^; otherwise, the convergence cannot be guaranteed 
[3l[6lll0]. However, for a problem with a complicated material interface, constructing a satisfactory body¬ 
fitting mesh is often costly, and this burden becomes more severe if the simulation involves a moving 
interface HU EH 125] because the mesh has to be generated repeatedly according to each interface location 
to be considered. In addition, some simulations, such as the particle-in-cell (PIC) method [211221120], can 
be carried out more efficiently on structured/Cartesian meshes. Due to these reasons, a wide variety of 
numerical methods based on Cartesian meshes have been developed. For an overview of these methods, 
we refer to [UlIilllOlllIllMllMl and the references therein. 

Immersed finite element (IFF) methods were recently introduced for solving interface problems. The 
main feature of IFF methods is that they can use meshes independent of the interface location, i.e., 
they allow interface to cut through the interior of elements in a mesh (see the last two plots in Figure 
[^. Hence, Cartesian (triangular or rectangular) meshes can be preferably employed in IFF methods 
to solve interface problems. We refer the readers to HU EH EH EEl EH for more features about IFF 
methods based on triangular meshes, and [17111111201122] for IFF methods based on rectangular meshes. 
We note that IFF methods in these literatures are applied in the continuous Galerkin formulation. 



Figure 1: from left: a simulation domain, a body-fitting triangular mesh, a non-body-fitting Cartesian mesh, and a 
non-body-fitting triangular mesh. 

The discontinuous Galerkin (DC) methods for elliptic boundary value problems can be traced back 
to 1970s (see mm) and they become increasingly popular recently as indicated by these survey articles 
and books [IIIIEIET]. Because there is no continuity imposed on the approximating function across 
the element boundary, DC methods can locally perform h-, p-, and hp- refinement flexibly and efficiently. 

For elliptic and parabolic equations, the interior penalty DC (IPDG) methods [H [T211221 SI] are well 
understood and widely used. The main feature of IPDG methods is that penalty terms are added on 
interior edges to stabilize the bilinear form of the scheme, so that the linear system is positive definite. 

In [221 EH], the IFF and IPDG ideas were combined together for solving interface problems on Cartesian 
meshes with local rehnement capability. To alleviate the issue of higher degrees of freedom in usual 
DC formulation, authors in m considered the so-called selective DG-IFF methods that employ DG 
formulation in selected elements while using the usual Galerkin formulation in the rest of the solution 
domain. Numerical examples have demonstrated that these DG-IFF methods can converge optimally, 
and our goal in this article is to theoretically establish the optimal a priori error estimates for DG-IFF 
methods that were discussed in ElEHlIHn]. 

The rest of the paper is organized as follows. In Section 2, we recall the DG-IFF methods originally 
proposed in |16[ll8j . In Section 3, we present a priori error estimates for these DG-IFF methods. An error 
estimate in a mesh-dependant energy norm is derived, and this error estimate is optimal according to the 
polynomials used in the IFF spaces. In Section 4, numerical experiments are provided to demonstrate 
features of DG-IFF methods. Brief conclusions are given in Section 5. 
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2 Discontinuous Galerkin Immersed Finite Element Methods 


In this paper, we adopt notations and norms of usual Sobolev spaces. For r > 1 and any subset G C fl 
that is cut through by F, we use the following function spaces: 

H^{G) = {vGH\G)-. ulonn^ GFF'(GnfI*),s = + or-}, H^{G) = H^{G) n H^{G), 

equipped with the norm 


From now on, we use G with or without subscripts to denote generic positive constants, possibly different 
at different occurrences, but they are independent of the mesh size and interface. 

Let {Th} with 0 < h < 1 be a family of triangular or rectangular Cartesian meshes of fl. An element 
cut through by the interface is called an interface element; otherwise, it is called a non-interface element. 
For each mesh Th, we let 7^® be the set of interface elements of Th and 7^" be the set of non-interface 
elements. We denote by £h the set of edges of Th- Also, let £h and £\ be the set of interior edges and 
boundary edges of Th, respectively. Similarly, if an edge is cut through by the interface, it is called an 
interface edge; otherwise, it is called a non-interface edge. Let and be the set of interface edges 
and non-interface edges, respectively. Moreover, we use and to denote the set of interior interface 
edges and interior non-interface edges, respectively. Without loss of generality, we assume that elements 
in Th satisfy the following conditions: 

(HI) If one edge of an element meets F at more than one point, then this edge is part of F. 

(H2) If F meets the boundary of an element at two points, then these two points must be on different 
edges of this element. 

For every interface element G 7^®, we assume its boundary intersects with the interface F at 
points D and E. Then, the line segment DE divides K into two sub-elements and K~ with 
K = U K~ U DE. With a given mesh Th on 17, we define the following broken Sobolev spaces: 

H\Th) = {v& L\n) : ViL G 7^,^!^ G H^K); 

ViLG7;\u|KGiLi(iL), v \ k ^ = +,-}. 

and 

HliTh) = {u G H^ijh) : = 0}. 

We now recall some standard notations for describing IPDG methods [gEain!. For each edge B, 
we associate a unit normal vector n^. If S G £h, we let Kb^i and Kb .2 be two elements that share B 
as the common edge and let be the outward normal with respect to 77^,1 • If i? G is taken 

to be the unit outward vector normal to d£l. For a function u defined on Kb.i U Kb. 2 , we denote its 
average and jump over i? G by 

{{■“S-B = 2{{'^\kb,i)\b + Mb = (u\kb,i)\b - {.u\kb,2)\b- 

If i? is a boundary edge, we set 

= Mb = w|b- 

For simplicity, we usually drop the subscript B from these notations if there is no danger to cause any 
confusion. 

To obtain a variational form for the interface problem 
function v G H^iTh), integrate both sides on each element K € Th, and apply the Green’s formula to 
have 

f fiTu-TvdX — f l3Tu-nKvds= f fvdX. ( 2 . 1 ) 

JK JdK JK 


(1.1) - (1.5), we multiply (1.1) by a test 
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Note that (2.11 holds regardless whether is a non-interface element or an interface element. For 
K G 7^", the derivation follows from the standard procedure. When K is an interface element, (2.11 
follows from applying the Green’s formula piecewisely over sub-elements of K determined according to 
the smoothness of u and v, then summing up over K and applying the flux continuity (1.5). Summarizing 
(2.1) over all elements we obtain 


f l3\7u-\7vdX— j • nsj-|z)] ds = j fvdX. 


( 2 . 2 ) 


^671 “ " B&Sh 

Since the solution u is continuous almost everywhere in fl, we can assume 

,.0 


^ H M ds = 0, Y, N 1^1 

B^Sh B^Sfi 


mI bl ds = 0, 


(2.3) 


for any constants e, a > 0, and > 0. Here \B\ denotes the length of B. Adding the two terms in 
(2.3) to (2.2), we obtain the weak form of interface problem (O - (|1.5|): Find u € such that 


ae{u,v) = {f,v), VvGHoiTh), 
where the bilinear form ae{-, •)■ Hh{^) x Hh{^) —>■ M is 


(2.4) 


{w,v) = 'Y [ /3Vw • VwdA — 'Y [ l,dVw • nsj-|t>] ds 


B^Sfi 


+<^Y i {{/3Vi;-nBS-Hds-k Y / M H ds, (2.5) 

and The weak form derived here for the interface problem ([^-([^ is in 

the same format as the standard weak form used in DG finite element methods for the usual elliptic 
boundary value problems miiiisi]. As suggested by DG finite element methods, the parameter e 
is usually chosen as —1, 0, or 1. Note that the bilinear form is symmetric if e = —1 and is 

nonsymmetric otherwise. 

We now introduce the IFF approximation of the broken space ido(7)i). For every element K G Th, 
denote by A^, i = 1, • • • , dif, the vertices of K. Here di^ = 3 or d/f = 4 depending on whether AT is a 
triangular or rectangular element. On each non-interface element K G Tj^, we let ipi, i = 1, - ■ ■ ,dx be 
the standard linear or bilinear hnite element nodal basis associated with the vertex Ai of K. The local 
FE space on K G is the defined as 


or bilinear [miSO] IFF 


Sh{K) = spanlipi : 1 < i < d^f}. 

On an interface element AT G 7),*, we let ())i, i = 1, • • • , dx be the linear 
nodal basis associated with vertex Ai. We let local IFF space on AT S 77( be 

Sh(,K) = : 1 < * < dx}- 

Then, we define the discontinuous IFF space over the mesh Th as follows: 

Sh{Th) = {VG : v\x G Sh{K)}, ShiTh) = {VG ShiTh) : v\an = 0}. 

One can easily see that Sh{Th) is a subspace of id?,(0). 

Finally, we state the DG-IFE methods for the interface problem ([^ - ([^ as: Find Uh G ShiTh) 
such that 

aeiuh,Vh) = if,Vh), yvhG ShiTh)- (2.6) 
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3 A Priori Error Estimation 


In this section, we derive a priori error estimates for the DG-IFE methods (2.61 in an energy norm 
II • ||/i : —>• K defined as follows 



1/2 


pVvVvdX+ 


|B|“ 


y] [w] ds 


(3.1) 


BgSh 


We first present a few lemmas required in the error analysis. By the standard scaling argument, one can 
show the following trace inequalities [37]: 

Lemma 3.1. (Standard trace inequalities) Let K he a triangle or rectangle with diameter h^, and B 
be an edge of K. There exists a constant C such that 

Mlhb) < C\Bp^\K\-^/^i\\vh2iK)+hK\\^v\\mK)), VvGHHK), (3.2) 

||Vi;|U2(b) < C\Bp^\K\-^^H\\^y\\LHK)+hK\PML-{K)), G H^(K), (3.3) 

where is the Hessian ofv. 

On an interface element K G Tf, we recall from UnilSSI that the local IFE space Sh{K) c H^{K). 
This implies that the trace inequality (3.2) is valid for v G Sh{K). However, since Sh{K) H‘^{K) for 
K&n in general, the second inequality ( |3.3| ) cannot be applied to functions in Sh{K). Nevertheless, in 
[32l|42], this trace inequality has been extended to IFE functions. We recall this result in the following 
lemma. 


Lemma 3.2. (Trace inequalities for IFE functions) Let Th he a Cartesian triangular or rectangular 
mesh and let K G Th be an interface triangle or rectangle with diameter hx and let B he an edge of K. 
There exists a constant C, independent of interface location but depending on the jump of the coefficient 
(3, such that for every linear or bilinear IFE function v defined on K, the following inequality hold 

||/3Vu • nB||i2(B) < O/i^^^^ll v^Vu||i2(^). (3.4) 

We now describe the interpolation with IFE functions. Eor K G Tff, the local interpolation operator 
is defined as IJ(x ■ ^'i^) Sh{K): 


d-K 

IlKn{X) = Yu[A{)f,,{X), KGTff. 

Eor K GTf, the local interpolation operator is defined as ^ : C{K) —>• Sh{K)-. 

dK 

/^,^u(X)=^u(H,)<)>.(X), KgT^. 

i=l 


On each non-interface element, we have the standard approximation theory for the finite element inter¬ 
polation: 

IW ~ Ih,Ky\\L^(K) + hxlu — Ih^Ky\H^{K) Si 'iKGTh- (3.5) 

On each interface element, the approximation property of the lEE interpolation proved in mm 
provides similar error bounds as follows: 


Ik ~ P,Ky\\L^{K) + hxlu — Ih,Ky\H^{K) < C!hx\\u\\ ^ 


(3.6) 


where the constant C is independent of interface location. For u G 
the interpolation defined by 


{IhU)\K = 


IlxU^ KGTff, 

^ ^ PI- 


let 4 : ^2(0) ^ Sh(Th) be 
(3.7) 
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Multiplying on both sides of (3.5) and (3.6), then summing up for all non-interface and interface 
elements, we can obtain an interpolation error bound on the domain Q as stated in the next lemma. 


Lemma 3.3. For u G satisfying the interface jump conditions (1.4) and (1.5), there exists a 

constant C such that 

( W'^ ~ ~ < C'/i||m||^2(q), (3-8) 

K&Th Ken 

where h = max hK- 
Ken 

The following lemma provides the approximation property of IhU in the energy norm || • ||?i. 

Lemma 3.4. Assume a < 1 in the energy norm (3.1). For euery u G F[^(fX), satisfying the interface 
jump conditions (1.4) and (1.5), there exists a constant C such that 

\\u-Ihu\\h<Ch\\u\\p[ 2 f^ay ( 3 - 9 ) 

Proof. By the definition of || • ||?,, we have 

\\u-hu\\l^ Y, f l3\Viu-hu)\^dX+ Y ^Wb-hujWh^sy (3-10) 

Ken-^^ Ben' ' 


For the first term on the right hand side, we use the estimate ( |3.8[ ) to have 

Y j hu)\^dX < Idmax Y l^maxh 


Ken 


Ken 




(3.11) 


where Pmax = max{/3 ,/3+}. Now we bound the second term in (3.10). Using the standard trace 
equality (3.2) and the approximation properties (3.5) or (3.6), we have 


\B\ 


|m IhU\ 11^2(5) < 

< 


~ IhU)\KB,i IIL2(b) + ||(m — IhU)\KB,2\\n(B)) 
ch-K \7 (ll^ - ^h^WlHKBn + II 

+ (ll^ “ ^/l^^lli2(ifB_2) + - IhU)\\\2(KB,2^ 

^ ^^KbJ\'^\\v{Kb,i) + ^^?ff3'^2ll^llnifi3,2) 

< Ch^ “ (ll^llv(iCi3,i) + Il^llv(ifi3,2)) ’ 

where V{K) = Fl^{K) for K G Tff and V{K) = H^{K) for K G Tf, and h = max hn- Also, the 


Ken 


second inequality is due to the shape-regular property of Cartesian triangular or rectangular meshes 
hKB,i ^ C'I^Ij * = 1,2. Thus, for a < 1, we get 


E ^III«-^HIIi^(B)<Ch 2 |kll| 2 (n)- 


B^Sh 


Finally, combining (3.11) and (3.12), we get (3.9). 


(3.12) 


□ 


The coercivity of the bilinear form ae(-, •) is analyzed in the following lemma. 

Lemma 3.5. There exists a constant k > 0 such that 

ae(vh,Vh) > KWvhWi, yvhGSh{Th) (3.13) 

holds for e = 1 unconditionally and holds for e = 0 or — 1 under the conditions that the penalty parameter 

7° 

’b 


CTg is large enough and a > 1. 
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Proof. From the definition of ae(-, •), we have 

a^ivh,Vh)= f /3\VvhfdX + {e-l) Y f • nsS K1 ds + Y i bhf ds. (3.14) 


B^Sh 


B^Sh 


We first note that, when e = 1, the coercivity follows directly from (3.14) and the definition of || • |[/j. 
If e = 0 or — 1, we need to bound the second term on the right hand side of (3.14). For each B G Sh, 
recall that Kb^i G Th, i = 1,2 are two elements sharing B as their common edge. If KB,i, z = 1 or 2 is 
a non-interface element, by the trace inequality (3.3) and inverse inequalities, we have 

\\{PVVh ■ nB)\KB,i\\L^{B) < /3max||(VUft,)|ifB 

^ ^Pma.xhj^^.W'S/VhWL^^KB.i) 

^Vh\\L^{KB,i)> 


/3n 


< c^^hA.\\^/pv^ 




‘-KB.i 


(3.15) 


where ^Smin = inin{/3“,/S+j, and /3max = max{/3“,/3+}. Then, using the assumption that a > 1 and 
by either the estimate ( 3.15[ ) or IFF trace inequality (3.4) depending on whether the element is a non¬ 
interface element or an interface element, we have 

f • tislluh] ds < |j • nB|||i2(B)|||uh] ||l2(b) 


< 


< 


1 
2 

C 1,-1 


(ll(/3Vu,, • nB)|ifB,ilU2(B) + ll(/3Vu/, • nB)\KB,2\\L-^(B)) II \\l^(b) 


< C 


1 


(lIv^Vu.ll 


2 


\\YP^Vh\\l2(^KB,2) 

Summing over all interior edges and using the Young’s inequality, we have 
Y f ■ nsl |uh] ds 

' Y (ll+ II^ 

B^Sh 




[w?l] ||l 2(B) 
[f/tl ||l2(B). 


B^Sh 


< 


^ ^ f E ^ 11 M ( Y (ll^/^vu,|| 


^|„/2ll h^B) 
IIV^Vu„|| 


1/2 


2 

L^{Kb.i 


LHKb,2), 


< 


E ^ E i^iiw 

B^Sh 


KeTh 

Then, for e = 0, we can choose 

and for e = —1, we can choose 


iL2(B)- 


(3.16) 


(5=1 and CTr > 


C 


(5 = - and a% > 2C. 


Substituting these parameters in (3.16) and then putting it in (3.14), we obtain (3.13). 


□ 


We also need an error bound for the IFF interpolation I^u on interface edges which has been proved 
in [32] • We present the result in the following lemma. 
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Lemma 3.6. For every u G H^{n), satisfying the interface jump conditions (1.4) and (1.5), there exists 
a constant C independent of interface location such that 


||/3(V(m - hu))\K ■ nsllisfs) < C 


(3.17) 


where K is an interface element and B is one of its interface edge. 

The assumptions of a in Lemma 3.4 and Lemma 3.5 suggest that we should choose a = 1 in our DG 
formulation (2.6). Now we are ready to prove the a priori error estimate for DG-IFE method (2.6). 

Theorem 3.1. Let u G he the exact solution to the interface problem <ini) to ( |1.5[ ) and Uh G 

Sh{Th) be the solution to (|2.6[) with a = 1, e = —1, 0, or 1. Then there exists a constant C such that 


\u-Uh\\h < ChWuWfjs^^y 


Proof. Subtracting the weak form (2.4) from the DG-IFE scheme ( |2.6[ ), we get 

ae{u-Uh,Vh) = 0, \fvhGShiTh). 


(3.18) 


(3.19) 


For every Wh G Sh{Th), using (3.19) and the coercivity (3.13), we have 
K\\uh-Wh\\l < ae{uh-Wh,Uh-Wh) =ae{u-Wh,Uh-Wh) 


< 


/ l3V{u-Wh) Wh)dX + i 113V{u-Wh)-nB^luh-Wh\ds 

Jk „ Jb 


Ken 


BgSh 


y^ / {{/3V {uh - Wh) • nsllu - Whj ds 

^ JB 


B^Sh 


= Ti+r2+T3+T4. 


E 

B^i-h 


\B[ 


(u - w/j] luh - Whj ds 


(3.20) 


We proceed to bound the terms Ti,i = 1, 2,3,4 in (3.20 1 . By the Cauchy-Schwarz inequality and Young’s 
inequality with parameter 5 > 0 , we can easily bound Ti and T 2 : 


1/2 


1/2 


Ti < 


E llv^V(n-n;,,)|| h(K) E \\y^V{uh - Wh)\\l2^K) 


\Ken 


\KeTh 


< 


1 


and 


.^/3max|lV(M 'U^/i)|lL2('f2) + ^ El {uh 'K!h)\\B2(K) 

Ken 

< C{S)\\V{u - Wh)\\l 2 ^Q) + 6\\uh - WhWl, 

T 2 < C(5) E + E 


(3.21) 


B^Sh 


B^Sh 


\B 


< C{5) E '-^\\miu-Wh)-nB}\\l2(B)+S\\uh-Wh\\l 


(3.22) 


B^Sh 


where C{5) emphasizes that this is a constant depending on 5. For T 3 , by the Cauchy-Schwarz inequality 
we have 

73 < E II - w?,) • nsj-|| l2(b)|| |m - w/i] ||l2(b). (3.23) 

B^Sh 
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First, using the standard trace equality ( |3.2[ ), we have 

II [u - Wft] ||i2(B) < \\{u-Wh)\KB,i\\L^(B) + \\{u-Wh)\KB,2\\L'^{B) 

^ {\W - '>^h\\L^(KB,i) +hKB.A\^^^-'^h)\\L^{KB,i)) 

+^^~KbI (II“ “ Wh\\L^(^KB,2) + ^ifB.2ll^(^ “ ^h)\\L^{KB,2)) ■ 


Then, by the trace inequalities (3.3) or (^3.4| , we have 


{{/3V(u/j - Wh) ■ nsS ||l2(b) < C (Vg'^^||\/^V(u/i - Wh)\\L^(KB,i) + {uh - Wh)\\L^(KB,2)) ■ 


Substituting the above two bounds into (3.23) and applying Young’s inequality, we obtain 


T3<C7(5)| ^ hj^Wu-WhWhi^K)^ ^ l|V(w-w;i)||i2(^) j +(5||Mft - w,,||^ 

Vi^GTh iCGTh / 


(3.24) 


For T 4 , we use the assumption a = 1, the Cauchy-Schwarz inequality, and Young’s inequality to have 

rO /• ^0 /• 


T, < Y. 


B^Sh 


^ ^ {u- w/j] lu-Whjds + s'^ I (Uh - [Uh - Wft,] ds 


4S\B\ 


IB 


\B\ 


IB 


(3.25) 


Again, by trace inequality (3.2), we have 

,.0 


|S| 


/* U 

jj,u-Whfds < (||(u-u;,,)|if^ J|^2(B) + ||(u-u>,,)|if^_J|^2(B)) 

- ^^~Kb,i ( 11 ^ “ ■^/i|Il 2 (Kb,i) + ^iCB,ilV(^ “ Wh)\\L'^(KB,i)) 

'^^^~Kb 2 {V-Wh\\L-^(KB,2) +hKB,2\\^{u-Wh)\\L-^{KB,2)f- (3-26) 


Using (3.26) in (3.25), we have 


T4<C{5)( hj/\\u-Wh\\l2i^K)+ Y \\'^(^-^h)\\l2(^K)] +Huh- 

KkgTh KeTh / 


WhWl- 


(3.27) 


Substituting (3.21), (3.22), (3.24) and (3.27) into (3.20) and choosing <5 = k/ 8 , we obtain 


\uh-Wh\\l < U ^ ||V(u-'u;?,)|||2(^)+C' ^ ^||{{/3V(M-u;,j)-nBS-|l|2(B) 


KeTh 

+C 'Y^ h-K \W ~ '^h\\L2f^j^-^ 

KeTh 


BGSh 


(3.28) 


linear or bilinear DG-IFE spaces (3.8) to have 


Now, we let Wh be the IFF interpolation I^u in (3.28) and use the optimal approximation capability of 


\un-hu\\l<Ch^u\\%^^^+Ch Y Y W'^{u-hu)-nB)\KBAh(By (3-29) 

BG£h *=1.2 


We now bound the second term on the right hand side of (3.29). If ATsy, i = 1 or 2 is a non-interface 
element, we use the trace inequality (3.3) to obtain 


II (/3V(u - hu) ■ nB)\KB,.\\Y b) < IIV(^^ - 4«)lli2(K,..) + hKB,. II V^u||i 2 (^^_,)) 

< C'hifg J|u||^2(^^ (3.30) 
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If Ks.i is an interface element, we use (3.17) to get 


\\{pv{u - hu) ■ nB)\KBA'L^B) ^ ■ (3.31) 

Due to the shape regularity of mesh Th^ we have the following bound on the union of interface elements 


Ker, 






^ hif < Ch\\i 


'H3{Q)- 


(3.32) 


iCGT' 


Summing up the estimates (3.30) and (3.31) over all interior edges, and using the bound in (3.32), we 
obtain 


E E ll(/3V(w - hu) ■ nB)|ifB,Jli2(B) < C'^l|w|l|3(n) + C'^lkll|2(a) < C'^l|w|l|3(n)- 

Be4 *=b2 


Then, substituting (3.33) to (3.29) we obtain 

\\uh — Ihu\\h < C!h\\u\\ 


Finally, the error estimate (3.18) follows from triangle inequality, (3.34) and (3.9). 


(3.33) 

(3.34) 

□ 


Remark 3.1. The DG-IFE methods proposed in this article and their related error estimation can he 
extended to arbitrary shape-regular unstructured interface independent meshes. 

Remark 3.2. The proof of Theorem 3.1 requires that the solution is piecewise H^, which is higher than 
the usual piecewise assumption imposed on methods using a finite element space based on linear 
polynomials. Consequently, our error estimate here is optimal according to the rate of convergence 
expected from linear polynomials but not with respect to the regularity of solution space. 


4 Numerical Examples 

In this section, we present numerical examples to demonstrate features of interior penalty DG-IFE 
methods for elliptic interface problems. Let the solution domain D be the open rectangle (—1,1) x(-l,l) 
and let the interface F be the ellipse centered at {xo,yo) = (— 0 . 2 , 0 .!) with semi-axes a = g^, b = |a. 
The interface separates 12 into two sub-domains, denoted by 12“ and 12“'', i.e., 

12 “ = {(a;, y) : r{x, y) < 1}, and 12+ = {(a;, y) : r{x, y) > 1}, 


where 




The exact solution u to the interface problem is chosen as follows 


u{x,y) = 


a'^b'^f^, ii{x,y)en , 

^b'^ -b ^ , if (a;, y) S 12 +. 


a 


(4.1) 


Here p is a parameter and we choose p = 5 in Examples 1-3 representing a solution with enough 
regularity, and p = 0.5 in Example 4 representing a solution with low regularity. Note that this solution 


does not satisfy the homogeneous boundary condition (1.2). We use this function for numerical verifi¬ 


cation because both the algorithm and the analysis in Section 2 and Section 3 can be extended to the 
nonhomogeneous boundary condition case via a standard treatment. 


Example 1: In this example, we present a group of numerical results for demonstrating the con¬ 
vergence of the DG-IEE methods on Cartesian triangular meshes. Additional numerical results on 
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rectangular meshes are provided in [161 UHl |20] . Specifically, the Cartesian triangular meshes {7^, h > 0} 
are formed by first partitioning fl into N x N congruent squares of size h = 2/N and then dividing each 
rectangle into two congruent triangles with one of its diagonal lines. 

First, we consider the case in which (/3“,/3^) = (1,10) representing a moderate discontinuity in 
the diffusion coefficient across the interface. The symmetric DG-IFE scheme is employed to solve the 
interface problem with parameters a = 1 and = 1000 for all interior edges. Errors of numerical 
solutions in the L°°, L^, and semi-H^ norms are reported in Table [ij Eor comparison, we also solve the 
same interface problem using the continuous Galerkin linear lEE method [551 HH] on the same meshes, 
and the corresponding numerical results are listed in Table [^ 


N 

■ 

rate 

• \\l2 

rate 

• \m 

rate 

10 

2.8333E-2 


3.7991E-2 


6.7917E-1 


20 

8.4503E;-3 

1.7454 

9.3605E-3 

2.0210 

3.4653E-1 

0.9708 

40 

2.5075E;-3 

1.7527 

2.3062E-3 

2.0210 

1.7456E-1 

0.9893 

80 

7.2318E-4 

1.7938 

5.6970£’-4 

2.0173 

8.7630E-2 

0.9942 

160 

2.0134E-4 

1.8447 

1.4140£’-4 

2.0105 

4.3903E-2 

0.9971 

320 

5.4720E-5 

1.8795 

3.5178E-5 

2.0070 

2.1972E-2 

0.9986 

640 

1.4450E-5 

1.9210 

8.7729E-6 

2.0035 

1.0991E-2 

0.9994 

1280 

3.7496E-6 

1.9463 

2.1903E-6 

2.0019 

5.4965E-3 

0.9997 


Table 1: Errors of linear DG-IFE solutions with /3 = 1, I3~^ = 10, = 1000. 


N 

■ 

rate 

• Ll2 

rate 

• \m 

rate 

10 

2.8619E-2 


4.4051E-2 


6.7876E-1 


20 

1.1416E-2 

1.3259 

1.1333E-2 

1.9587 

3.4808E-1 

0.9635 

40 

5.3027E;-3 

1.1062 

2.8882E-3 

1.9722 

1.7641E-1 

0.9805 

80 

1.9396E-3 

1.4510 

7.3078E-4 

1.9827 

8.9155E-2 

0.9845 

160 

1.0689E-3 

0.8596 

1.8726E-4 

1.9644 

4.5387E-2 

0.9740 

320 

5.4774E-4 

0.9646 

5.1220E-5 

1.8702 

2.3305E-2 

0.9617 

640 

2.7498E-4 

0.9942 

1.6771E-5 

1.6107 

1.2277E-2 

0.9247 

1280 

1.4113E-4 

0.9623 

7.1759E-6 

1.2248 

6.7321E-3 

0.8668 


Table 2: Errors of the classic Galerkin IFE solutions with (5 =1, Z?"’' = 10. 


The data in Table [^ clearly demonstrate that the convergence rate of the DG-IFE method in the 
semi-H^ norm is optimal which corroborates the a priori error estimates (3.18) for DG-IEE methods 


since the semi-H^ norm is part of the energy norm. In addition, the data in this table indicate that 
the convergence rate of the DG-IEE method in the norm is also optimal. However, from Table 
we can see that the convergence rates of the Galerkin IFE solution in the and norms start to 
deteriorate when the mesh size becomes smaller than h = 2/320. This comparison indicates that the 
DG-IEE methods are more stable than the continuous Galerkin lEE method. 

Next we consider the case involving a larger discontinuity in the diffusion coefficient, z.e., /3~ = 1, 
and /?■*■ = 1000. We use nonsymmetric DG-IFE scheme for this experiment and choose erg = 1000 for 
all interior edges. As demonstrated by the data in Table the DG-IEE solutions converge optimally in 
the and semi-71^ norms. 


Example 2: From Tables[2and[^ it is interesting to note that errors in the DG-IFE solutions gauged 
in norm are much smaller than those in the classic Galerkin IFE solutions when the mesh size is 
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Pointwise Error 


N 

■ 

rate 

• L2 

rate 

• \m 

rate 

10 

1.9381F-2 


1.5338F-2 


2.1012E-1 


20 

1.2420F-2 

0.6420 

6.0704F-3 

1.3373 

1.3141F-1 

0.6772 

40 

4.0332F-3 

1.6227 

1.4957F-3 

2.0210 

6.9522F-2 

0.9185 

80 

9.9934F-4 

2.0129 

3.6124F-4 

2.0498 

3.5490F-2 

0.9701 

160 

2.7965F-4 

1.8374 

8.9863F-5 

2.0072 

1.7949F-2 

0.9835 

320 

8.0700F-5 

1.7930 

2.1864£’-5 

2.0392 

9.0223F-3 

0.9924 

640 

2.2017F-5 

1.8740 

5.3914F-6 

2.0198 

4.5229F-3 

0.9963 

1280 

5.9615F-6 

1.8848 

1.3343F-6 

2.0146 

2.2641F-3 

0.9983 


Table 3: Errors of linear DG-IFE solutions with j3 =!,/?+ = 1000, = 1000. 

small enough. It has been observed that the classic IFE solution has a so called “crown” shortcoming as 
demonstrated by the plot on the left side of Figure meaning its point-wise accuracy is much worse in 
the vicinity of the interface than the rest of the solution domain. We think this severity of inaccuracy is 
caused by the discontinuity in the IFE functions across the interface edges. Nevertheless, the DG-IFE 
methods contain penalty terms that can alleviate the adverse impacts from the discontinuity across ele¬ 
ment edges, especially those from interface edges. Therefore, DG-IFE methods can usually outperform 
the classic Galerkin IFE method around the interface as demonstrated by the plot on the right side of 
Eigurej^ lEE solutions in these plots are generated on a mesh formed by partitioning il into 160 x 160 
congruent rectangles first, then generating triangular elements by the diagonal line of these rectangles. 



Figure 2: Point-wise error of Galerkin IFE solution and DG-IFE solution. 

Example 3: One desirable feature of the DG formulation is the local adaptivity in mesh or polynomi¬ 
als because this formulation does not require the inter-element continuity of finite element functions. The 
combination of IFE spaces and the DG formulation leads to a new class of finite element methods that 
allow local mesh refinement while maintaining the possibility of using the desirable structured Gartesian 
meshes for solving problems with nontrivial interface geometry. This example is for demonstrating this 
feature of the DG-IFE methods. 

The discontinuity in the coefficient of the interface problem limits the smoothness of the exact 
solution around interface. The low regularity of the solution around interface usually has an adverse 


12 






















impact on the accuracy of the numerical solution around the interface. To overcome this challenge, one 
can employ more finite element functions around the interface and a way to achieve this is to refine 
interface elements. To demonstrate this idea, we consider the same example described above for a larger 
coefficient jump (/3“ = 1, = 1000). We start with a uniform Cartesian mesh consisting of 

10 X 10 rectangles. For A: > 1, the mesh is generated by refining the previous mesh via 

cutting each of its interface elements into four congruent rectangles by connecting midpoints on opposite 
edges. We then solve the interface problem 0) 


(1.5) on a sequence of 6 such meshes generated by 
the refinement procedure described above using the nonsymmetric DG-IFE method with tr^ = 1000 on 
all internal edges. Errors in the and semi-iJ^ norms generated on each mesh are presented in 


Table In the second column, |7)| 


Wi 


denotes the number of elements in the mesh 7)) 


(fc) 


The numbers 


of degrees of freedom are listed in the third column. The initial mesh 7), 


(0) 


and Tf, 


( 4 ) 


are illustrated in Figure]^ 


and the refined meshes 7)^^^^ 





Figure 3: Solution meshes . 


Mesh 

IT-Wl 

1 >h 1 

DoF 

II ■ \\ l °° 

• \\d 

• HI 

T-(0) 

'h 

100 

400 

1.5064E-2 

1.5667E-2 

1.9393E-1 

'h 

178 

712 

1.9242E-2 

8.7358E-3 

1.5041E-1 

-7"(2) 

'h 

334 

1336 

1.2799E-2 

5.9608E-3 

1.2761E-1 

-7"(3) 

'h 

646 

2584 

1.2789E-2 

5.7452E-3 

1.2321E-1 

-7"(4) 

'h 

1258 

5032 

1.2484E-2 

5.6245E-3 

1.2233E-1 

q~{^) 

'h 

2470 

9880 

1.2431E-2 

5.5988E-3 

1.2213E-1 

'7"(6) 

>h 

4882 

19528 

1.2424E-2 

5.5940E-3 

1.2209E-1 


Table 4: Errors of NIPDG-IFE solutions on meshes with local rehnement 

The data in Tablej^show that the global error in the and semi-7f^ norms are significantly reduced 
in the first two steps of local mesh refinements. But further refinements performed on interface elements 
do not reduce the global error as much as the first two steps. We believe this is because, after the first 
two refinements, errors of the DG-IFE solutions over non-interface elements become more significant. To 
further increase the accuracy of DG-IEE solutions, refinement on non-interface elements is also necessary. 
To demonstrate this idea, we simulate an adaptive local mesh refinement over the whole solution domain. 
Since a posteriori error estimators for IFE methods are not available yet, we use the actual error as the 
“ideal” error indicator to guide the local refinement just for a proof of concept. 

We start from a uniform Gartesian mesh 7)^^°^ consisting of 10 x 10 rectangles. Eor fc > 0, we produce 
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a DG-IFE solution Uh for the interface problem on the mesh 7^^^^ and compute the local semi-iJ^ norm 
error \u — Uh\H^{T) on each element of We sort these local errors from the largest to the smallest 

and use this order to form the smallest collection ^ of the first few elements such that 

|u — > 0 |m — Uft,|^i(Q). (4.2) 


Then, we generate a new mesh 7),^*^^^ by refining via cutting each of the elements in into four 
congruent rectangles by connecting midpoints on its opposite edges. The computation repeats over the 
new mesh ■ 

We choose 9 = 0.2 and conduct adaptive DG-IFE scheme on each locally refined mesh. In Figure 


1^ we show the initial mesh the refined meshes and and errors of the DG-IFE 

solutions generated on those meshes. From these plots, we can easily see that the global error in the 
DG-IFE solution is reduced as the local mesh refinement automatically deploys smaller elements at 
locations needed according to the “ideal” error indicator while maintaining the Gartesian structure of 
the meshes. 

To see the effectiveness of the adaptive DG-IFE methods, we compare the errors in DG-IFE solutions 
generated via adaptive mesh refinement to errors of DG-IFE solutions generated on uniform meshes with 
comparable degrees of freedom in Figurej^ The semi-71^ norm errors presented in the right plot in Figure 
shows that the magnitude of errors in adaptive DG-IFE method are smaller than in the method with 
a uniform mesh when their degrees of freedom are comparable. However, the order of convergence for 
both schemes are optimal by comparing their errors with the reference line of slope —1/2 (the same 
criteria is used in 0)- We note that the exact solution u defined in ( |4.1[ ) with p = 5 is piecewisely 
u G 77^(11) (see the left plot in Figure^, although the global regularity is impacted by 


smooth, i.e. 

the discontinuity of the coefficients. The optimal convergence of the numerical solution in uniform mesh 
refinement confirms our theoretical error estimate (3.18) in Section 3. 





Figure 4; Point-wise error of DG-IFE solutions on meshes: 


Example 4: In many applications, because of the insufficient regularity in the data, solutions to 
the involved boundary value problems may not be piecewisely smooth enough for a certain convergence 
theorem to hold. In such cases, finite element or DG methods based on uniform mesh refinement usually 
fail to converge optimally, but adaptive FE or DG methods with suitably designed mesh refinement 
strategies can still generate optimally convergent numerical solutions 0IH1. For interface problems with 
discontinuous coefficients, the phenomenon is similar. In this example, we demonstrate behaviors of 
DG-IFE method in adaptive and uniform mesh refinements for solving interface problems whose exact 
solution has singularity. In particular, this example indicates that the DG-IFE method with adaptive 
refinement on interface independent meshes can satisfactorily handle interface problems whose exact 
solutions are less smooth. 

We choose the exact solution u in the form of (4.1) with p = 0.5 such that it does not satisfy the 
regularity condition required by Theorem 3.1 on the convergence of the DG-IFE method. This lack of 
regularity is caused by the singularity of the exact solution at the center of the ellipse as depicted in the 
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DG-IFE 


10 ° 



10 ^ 10 ° 10 * 10 ^ 

# Degree of Freedom 


Figure 5: Left: the exact solution u in Example 3. Right: a log-log scale plot of errors of in numerical 
solutions generated by the adaptive and uniform mesh refinement. The red line with a —0.5 slope is for 
the reference of optimal convergence. 


left plot of Figure]^ The coefficients are set to be /3“ = 1, /3^ = 10. First, we present the data generated 
by the DG-IFE method with uniform mesh rehnement in Table Specihcally, these data are produced 
by the nonsymmetric DG-IFE method with the penalty = 100 on every edge. Errors in 
and semi-R^ norms are reported, and the convergence rates of the DG-IEE method are obviously not 
optimal in all three corresponding norms. 


N 

■ 

rate 

■ Ll2 

rate 

1 • \h^ 

rate 

10 

4.2599E-2 


1.8474E-2 


1.0415E-1 


20 

5.7301E-2 

-.4278 

8.0329E-3 

1.2015 

5.6787E-2 

0.8751 

40 

4.5270F;-2 

0.3400 

5.5008E-3 

0.5463 

4.2830E-2 

0.4069 

80 

3.5386F;-2 

0.3554 

3.8353E-3 

0.5203 

3.2011E-2 

0.4201 

160 

2.7412E-2 

0.3684 

2.6964E-3 

0.5083 

2.3774E-2 

0.4292 

320 

2.1075E-2 

0.3793 

1.9024E-3 

0.5032 

1.7573E-2 

0.4361 

640 

1.6098E-2 

0.3886 

1.3441E-3 

0.5012 

1.2940E-2 

0.4415 

1280 

1.2237E-2 

0.3967 

9.5014E-4 

0.5005 

9.5021E-3 

0.4461 


Table 5: Errors of bilinear nonsymmetric DG-IFE solutions with p = 0.5 


Next, we report the performance of the adaptive DG-IFE method for solving the same interface 
problem with this less smooth exact solution. As in Example 3, we use the exact error as an “ideal” 
error indicator for mesh refinement just for a proof of concept. Starting with a uniform mesh 
consisting of 10 x 10 rectangles, we perform the local mesh refinement based on the same rule as the 
one in (4.21 with the threshold 6 = 0.2. Errors in semi-R^ norm are depicted in Figure]^ in which, as a 
comparison, errors from uniform mesh refinement are also plotted. It is obvious that the adaptive DG- 
IFE method is far more accurate than the DG-IFE method based uniform meshes when their numbers 
of degrees of freedom are comparable. Moreover, comparing the errors with the reference line of slope 
— 1/2, it is obvious that the rate of convergence of the adaptive DG-IFE method is close to optimal 
from the point view of the degrees of freedom while the rate of convergence of DG-IFE method based 
on uniform mesh is not optimal. Some meshes in the process of the local refinement are presented in 
Figure from which one can observe that the refinement is around the center of the ellipse where the 
exact solution is singular. 
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# Degree of Freedom 


Figure 6: Left: the exact solution u in Example 4. Right: a log-log scale plot of errors of in numerical 
solutions generated by the adaptive and uniform mesh refinement. The red line with a —0.5 slope is for 
the reference of optimal convergence. 
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Figure 7: Solution meshes: and 


5 Conclusion 

In this article, we establish a priori error estimates for interior penalty discontinuous Galerkin methods 
with immersed finite element functions for elliptic interface problem. The method can be used on 
Cartesian meshes that are independent of the interface. The analysis here shows that the order of 
convergence of these DG-IFE methods is optimal in the energy norm from the point of the polynomial 
degree in the finite element spaces. With the enhanced stability, these DG-IFE methods outperform the 
the classic Galerkin IFE methods, especially in vicinity of the interface across which the exact solution is 
usually less smooth. The proposed DG-IFE schemes allow efficient local mesh refinement while preserving 
the Cartesian structure of meshes provided that a posteriori error estimators are available. 
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